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When a nematic liquid crystal is confined in a porous medium with strong anchoring conditions, 
topological defects, called disclinations, are stably formed with numerous possible configurations. 
Since the energy barriers between them are large enough, the system shows multistability. Our lat- 
tice Boltzmann simulations demonstrate dynamic couplings between the multistable defect pattern 
and the flow in a regular porous matrix. At sufficiently low flow speed, the topological defects are 
pinned at the quiescent positions. As the flow speed is increased, the defects show cyclic motions and 
nonlinear rheological properties, which depend on whether or not they are topologically constrained 
in the porous networks. In addition, we discovered that the defect pattern can be controlled by con- 
trolling the flow. Thus, the flow path is recorded in the porous channels owing to the multistability 
of the defect patterns. 



Topological defects are observed in various types of 
phases associated with spontaneous symmetry breaking 
such as crystals and ferromagnetic materials pQ. In 
the early stage of the phase transitions, many defects 
are created, and their total amount then decreases with 
time to reduce the free energy [2J. If frustrations are 
internally included or imposed by external fields, the 
defects are sometimes stabilized. On the other hand, 
re-organization of topological defects is often important 
under non-equilibrium conditions. For instance, pin- 
ning and depinning of quantized vortices in a type II- 
superconductor induce electrical resistance [3J. Motions 
of quantized vortices also lead to mutual frictional inter- 
actions in a flowing supcrfluid [2 . 

A nematic liquid crystal (NLC) is an ideal system for 
studying behaviors of topological defects because of its 
slow dynamics and large lengthscale |3H6]. When NLCs 
are confined in cavities [7J [5] , microfluidic devices [51 ITU] 
and porous media [11H17] . the anchoring interaction be- 
tween the director and the solid surface imposes frustra- 
tions on the elastic field. In a complex geometry, many 
topological defects in NLCs are thus sustained for a long 
time [H [T5H2T] . A NLC in a porous material exhibits 
slow glassy behavior and the resulting memory effects 
[151 1161 119j . In a previous paper, we studied the mem- 
ory effect by Monte Carlo simulations, focusing on the 
role of the bicontinuity of the matrix (2JJ . After cooling 
from an isotropic state, many disclination lines running 
through the channels of the porous network remain, with 
numerous possible trajectories. Because each trajectory 
is at one of the local energy minimum states and the 
energy barriers among them are often larger than the 
thermal energy, the system shows non-equilibrium be- 
haviors similar to those in spin glasses. We found that 
the memorization is attributed to the re-configurations 
of the defect structure. The applications of an electric or 
magnetic field can bring the system into a new different 
state, and the new configuration remains even after the 
field is removed. Here, irreducible disclinations, which 
are topologically concatenated to the solid matrix, tend 
to enhance the memory effect. 



In this Letter, we study a NLC flowing in a regular 
porous medium. In particular, we focus on dynamic cou- 
plings between the flow and the nonergodic disclination 
pattern. We hope our findings will facilitate research on 
microfluidics using liquid crystals [SJ [JO] . 

First, we prepare a pattern for a porous material by in- 
troducing a scalar variable in a cubic lattice according 
to Rcf. [2"T] . </> varies smoothly in space. Then, we parti- 
tion the cubic lattice into three portions by using 4>: V is 
the ensemble of lattice sites representing the solid mate- 
rial of the porous matrix; M denotes the nematic fraction 
and S is the ensemble of lattice sites at the interface. For 
the portions M and S , we introduce a nematic tensorial 
order parameter Q a p [22]. We employ the Landau— dc 
Gennes free energy density [23], 
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where op, /3p, 7f, L and Ae are phenomenological ma- 
terial constants. E is an external field. The repeated 
subscripts indicate summations over Cartesian indices. 
The free energy functional of the system is given by 

FiQap} = / drf LdG (Q a/3 ) - w / drQ a pd a dp,(2) 
JM+s Js 

where d is the normal vector of the porous matrix, 
d = V0/|V</>|. The second term represents the molec- 
ular anchoring interaction between the surface of the 
porous matrix and the nematic orientational order, and 
w is its material constant. We assume w > 0, i.e., that 
homeotropic anchoring is preferred. 

The hydrodynamic flow u is determined by a momen- 
tum equation [23j . 
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where p is the material density and p is the pressure. 
Here, a^o is the distortion stress tensor, which is given 
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Further, a' a p is the viscous stress tensor, which is given 
by 

a 'a/3 —PlQaf}Q^vA^ v + fi±A a p + feQapApp + (3 6 Q p^A^a 
+ 7^2^/3 - VlQa^N^p + HlQpftNfia, (5) 

where A a p = (S/ a up+\7 pu a )/2 is the symmetric velocity 
gradient tensor, and N a p is the time rate of change of 
Q a /3 with respect to the background flow. The latter is 
defined as 
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where ui a p = (V a u^ — V pu a )/2 is the asymmetric veloc- 
ity gradient. In Eq. (i = 1, 4, 5, 6) and /ij (i = 1, 2) 
are phenomenological parameters having the dimension 
of viscosity. The time development of the nematic order 
parameter is given by 
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The above nemato-hydrodynamic equations are solved 
with non-slip boundary conditions at the surfaces of the 
porous matrix. 

Lattice Boltzmann simulation is one of the most effi- 
cient numerical methods of investigating flow behaviors 
in a complex geometry [24]. The bounce-back method 
realizes non-slip boundary conditions on the velocity 
fields at the surfaces of the porous matrix [23]. To our 
knowledge, two lattice Boltzmann algorithms for simu- 
lating nemato-hydrodynamics have been proposed |25l - 
128] . In this study, we adopt a numerical scheme devel- 
oped by Care and coworkers. Here, we omit the details 
and follow their notations (see Refs. [251 123)- We per- 
formed three-dimensional simulations with four speeds 
and coordination number 27 (D3Q27). The weight fac- 
tors of the four speeds U are given by ti — 8/27 for 
d = (0,0,0), U = 2/27 for a = (±1,0,0), U = 1/54 for 
d = (±1,±1,0) and U = 1/216 for c t = (±1,±1,±1) 
with appropriate permutations. Here, c t is the non- 
dimensional velocity vector. 

In our simulations, we impose an external force density 
/ to flow the NLC. Then, we modify Eq. (31) in Ref. [26] 
to 
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where Y]j represents summation over 27 velocity vectors 
Cj. We employ a F = —0.133, /3p = 2.133, 7f = 1.6, w = 
0.1, L = 0.1, Ae = 1.0, m = 3.0, fi 2 = -1.0, t = 1.1, 



f]o = 0.3, and 772 = 0.5. This set of parameters gives a 
correlation length £ = 0.87A, an anchoring penetration 
length d a = A, and a bulk nematic order = 0.37, 
where A is the unit cell length. The resulting viscosity 
ratio of this anisotropic liquid is t]x/t)\\ — 1-1. Here 7711 
and rj± are the effective viscosities when the director is 
parallel to the flow and parallel to the velocity gradient, 
respectively. We adopt a bicontinuous pattern with cubic 
symmetry, which we have studied and is denoted as BC 
in Ref. [21]. The size of the unit cell of the bicontinuous 
cubic is L = 32 A. We use a simulation box of 64 3 with 
periodic boundary conditions. 

In NLCs, the flow speed is characterized well by the Er- 
icksen number Er, which is the ratio of the viscous and 
the elastic forces. It is given by Er — \ f\£ 3 /K, where K 
is the Frank elastic modulus expressed as K = 9LQ^/2 
in a one-constant approximation, and £ is the character- 
istic length of the distortion. Because it is difficult to 
determine I quantitatively in a complex geometry, we as- 
sume that I is of the order of the radius of the narrowest 
channels, namely, I = L/4. 

We quenched the system from an isotropic state to 
a nematic state without any external field E or bulk 
force /. Because of the anchoring effect, the director was 
strongly deformed, therefore, disclination lines remained 
even after long duration annealing. Figure [l|a) shows 
one of the stable defect configurations, where the defect 
positions are represented by the isosurfaces of the elastic 
energy density / c i = i(V M Q Q( g) 2 /2. Here, the remaining 
disclinations run through the channels rather randomly. 
They are closed without ends, and some are concatenated 
to the solid matrix. Then, we applied a strong external 
field pulse along the z-axis (E z = 0.1). After the field 
was removed, two types of disclination loops formed in 
the x-y planes, as shown in Fig.[ljb). Half of the remain- 
ing disclinations (red loops in Fig. [ijb)) are localized at 
the narrowest sections of the channels. The line defects 
locally have the topological charge s = —1/2, so the loops 
are reducible continuously into hyperbolic point defects 
with s = — 1. On the other hand, the other half of the 
loops (blue loops in Fig. [ITb)) consists of disclinations, 
which locally have the topological charge s = 1/2. How- 
ever, they are irreducible to point defects, because the 
solid necks of the BC penetrate the loops. Hereafter, the 
former and latter defect groups are referred to as type-I 
and II, respectively. This structure is the same as that 
obtained in the Monte Carlo simulations [5T] . It also ap- 
pears that this configuration is probably one of the global 
energy minimum states, which are also degenerated along 
the x and y-axes. We employ this configuration as an ini- 
tial state. 

First, we impose force f z to flow the NLC toward the 
z-direction. In NLCs, defect motion is caused by not only 
hydrodynamic convection, but also rotations of the direc- 
tor field. If the imposed force is weak enough, both types 
of defects show only small displacements, which are de- 
termined by the balance between hydrodynamic convec- 
tion and director rotation (not shown here). As the flow 
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speed is increased, the director rotation can no longer 
compensate for the hydrodynamic convection, so the de- 
fect pattern becomes unstable against the flow. Figure [2] 
shows the temporal changes in (a) the defect pattern and 
(b) the cross section of Q zx at an x — z plane that cuts 
the centers of the necks of the BC. The applied force 
strength is f z = 0.006, which corresponds to Er = 50. 
Under this force, the type-I defects start to exhibit cyclic 
behavior. As shown in Fig.pfb), the director field is grad- 
ually distorted with time. When the elastic distortion 
energy becomes comparable to the anchoring energy, the 
anchoring at the downstream side of the porous surface 
is transiently broken (see t — 270 in Fig. [2]jb)). When 
the anchoring breaks, new defects are born from the sur- 
face. Then, the new defects and the old type-I loops fuse 
into new type-I disclination loops. During this cycle, the 
type-II loops keep to be at the quiescent positions. With 
a further increase in the flow speed, the type-II loops 
also become unstable. Figure [2^c) shows snapshots of 
the defect pattern at f z — 0.013. Because they are con- 
catenated to the channels, they cannot move without a 
topological transformation. The type-II defects repeat 
another cycle, in which they are deformed, scissored and 
reconnected sequentially. 

In Fig. [3]ja), we plot the temporal changes in the av- 
eraged orientational order along the z-axis (Q zz ) . When 
the force is weak (f z < 0.004), the orientational order 
is almost fixed at the same value as that at rest. Sub- 
ject to an intermediate force (0.006 < f z < 0.008), the 
orientational order shows an oscillation mode, which cor- 
responds to the repeated cycle of the type-I defects shown 
in Figs. [2|a) and (b). For strong force (f z > 0.010), two 
oscillation modes are observed in (Q zz ). The mode of 
the shorter period is the same as that for f z = 0.006; 
on the other hand, that of the longer period corresponds 
to the cycle for the type-II defects shown in Fig. [2|c). 
The frequencies of the two oscillation modes are shown 
in Fig. ^b). Both frequencies decreased with decreas- 
ing Er. It is likely that they show critical behaviors as 
w a \Er — Er c \ k with k = 1.0. The critical values of Er 
are estimated to be Er\ = 26 for the type-I defects and 
Er l J = 64 for type-II defects by extrapolating straight 
lines. However, the true critical points are not accessi- 
ble because of our simulation's constraints. Therefore, 
we cannot yet conclude whether these bifurcations are 
continuous or not, and determine the exponent k. 

Figure [3tc) shows the temporal changes in the flow 
speeds (ill) averaged inside the porous matrix. Above 
the critical strength for the cyclic defect behavior, the 
flow speed exhibits oscillation, the period of which co- 
incides with that of (Q zz ). In Fig. [3](d), we plot the 
apparent viscosity rj z , which is obtained by averaging it 
over 1000 < t < 2000. Here, the viscosity is normalized 
by that of an isotropic liquid ?7i SO , which is determined by 
the same simulation method. Below the critical force, the 
apparent viscosity decreases marginally with Er. That is, 
the NLC exhibits a shear-thinning rheology. At the two 
critical values of Er for the onset of the oscillations, small 



kinks in the apparent viscosity are observed, indicating 
that the defect pattern influences the flow properties. 

Next, we impose an external force toward the in- 
direction. Under a weak force, all the defects are pinned 
at their original positions, showing marginal elongations 
along the flow (not shown here). Above a certain force, 
we found that the defect pattern switched to an orienta- 
tion along the flow direction. The switching process of 
the defect pattern under f x — 0.005, which corresponds 
to Er = 42, is shown in Fig. |4j Once the orientation is 
switched, the new defect pattern is sustained owing to its 
multist ability even after the flow is stopped. The appar- 
ent viscosity for the x-oriented flow, r/ x , is also plotted in 
Fig. [3^d) . Below the critical Ericksen number for switch- 
ing, the apparent viscosity along the x-axis is larger than 
X] z . At the critical Ericksen number (Er c j_ = 25), rj x 
decreases abruptly to r] z . After switching, the system ex- 
hibits the same oscillations as those for the parallel flow 
in Fig. § 

Note that the oscillation frequencies and the cor- 
responding Ericksen number depend on the anchoring 
strength, because the anchoring breaking plays a cru- 
cial role in the cyclic behaviors of the topological defects. 
Thus, the values displayed are not universal for any an- 
choring strength and pore size. Moreover, qualitatively 
distinct behaviors may be observed in extreme cases. For 
instance, defects may be transported unboundedly with- 
out the cyclic motions if the anchoring strength is in- 
finitely large, although our simulations cannot access this 
regime. Because the Ericksen number does not factor 
in the anchoring effect, it cannot describe all the flow 
properties. If the morphology of the porous matrix is 
the same, the director pattern is characterized by a non- 
dimensional parameter W = wL/K in the absence of 
flow. Two dimensional mapping of the non-equilibrium 
behaviors onto a W — Er plane will help to understand 
the flow properties of NLCs in confined systems. We need 
to study them more quantitatively. 

We studied the flow behaviors of a NLC in a porous 
medium by lattice Boltzmann simulations. In the case 
of strong anchoring, the pattern of the director field can 
adopt many (meta-)stable configurations. We found that 
this multistability leads to peculiar flow properties of the 
confined NLCs. The dynamic behaviors of topological 
defects under flow depend on whether they are topologi- 
cally concatenated to the matrix. 

We also found that the apparent viscosity depends on 
the mutual direction of the flow; i.e., the NLC flows eas- 
ily along the orientation of the averaged director field. 
Imposing an intermediate flow along an incommensurate 
direction switches the global director field to the new di- 
rection. Upon this switching, the apparent viscosity de- 
creases abruptly to that for the commensurate flow. Be- 
cause this defect pattern is maintained even after the flow 
is stopped, we can say the flow path can be recorded in 
the porous network. As reported in a previous paper, the 
emergence of topological defects enhance the efficiency of 
memorization [5T] . We also expect that the flow path in 
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a porous matrix can be dynamically selected by changing 
the multistable defect patterns. In our simulations, the 
difference between the apparent viscosities of the parallel 
and perpendicular flows is not large compared to that in 
Ref. |26j . In other words, our system is more isotropic 
in viscosity. This small viscosity difference is attributed 
to the employed parameters, which are restricted in our 
present scheme. When a NLC having a large viscosity 
anisotropy is used, the difference in the apparent viscos- 
ity will be more remarkable. 

We employed a bicontinuous matrix with cubic sym- 
metry as an example. Because the configuration of topo- 
logical defects depends on the morphology of the net- 
work, our results cannot be directly applied to any ma- 
trices. However, we consider that dynamic coupling be- 
tween the multiple stabilities and the flow will be found 
ubiquitously for NLCs confined in any porous network. 
We also hope our results will improve the understand- 
ings of the flow properties of NLCs in porous media and 
microfluidic channels. 
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FIG. 1. (a) Snapshot of disclination lines of nematic liquid 
crystals in a porous medium after zero field cooling from an 
isotropic state. Solid green objects represent the porous ma- 
trix, and red curves are the remaining disclinations. (b) De- 
fect pattern remaining after an external field pulse is applied 
along the z-axis. Red and blue curves represent reducible 
(type-I) and irreducible (type-II) disclinations, respectively. 
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FIG. 2. Time evolutions of (a) the defect pattern and (b) 
Q zx at an x-z plane under f z = 0.006. In (b), the black areas 
are the solid matrix, and white crosses indicate the defect 
positions. At t = 270, the anchoring condition is broken, (c) 
Time evolution of the defect pattern under f z = 0.013. The 
type-II (blue) defects also exhibit cyclic behavior. 
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FIG. 3. (a) Temporal changes in the averaged remnant order 
(Qzz). (b) Frequencies of the cyclic motions of the type-I and 
II defects as a function of the Ericksen number. Broken lines 
are provided for better visualization, (c) Temporal changes in 
the flow speed {u z ) of the liquid crystal in the porous matrix. 
The force density is increased from f z = 0.002 to 0.01. (d) 
Effective viscosities for the parallel and perpendicular flows 
plotted with respect to Er. The viscosity is normalized by 
that of an isotropic liquid ?7i so . At the arrows, the defect 
patterns are re-organized. 
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FIG. 4. Switching behaviors of topological defects under a 
perpendicular flow. A force f x = 0.005 is imposed for the 
defect pattern aligned along the z-axis. 



